clear
capture 
clear matrix
mat drop _all
set more off
set matsize 1000
*change here to relevant directory"
global tree "C:\Users\tbold\Dropbox\Tessa and Tobi\JEEA Replication"
************************************************************************************************************
*****this creates the data to be used for comparing fit of coalitional and individual deviations model******

/* Standards in the literature (Townsend, Ligon, Laczo)
1) use only 1976-1981 because of concerns about measurement of income in later years
2) use adult equivalent consumption
3) use only households that are there every period, i.e. a balanced panel
4) consider each village separately
*/
global data	"$tree\Data"
global latex "$tree\Latex"
global programmes "$tree\Programmes"

*global data "C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\Data"

global input 	"$data\Input"
global output 	"$data\Output"
global dofiles 	"$data\Dofiles"

cd 	"$data"

matrix MOMENTS=J(4,4,0) //1-3 Aurepalle, Kanzara, Shirapur, 4 ALL villages
matrix MOMENTSRAW=J(4,4,0) 
matrix MOMENTSDIFF=J(4,4,0)
matrix COV_MOMENTS_1=J(4,4,0)
matrix COV_MOMENTS_2=J(4,4,0)
matrix COV_MOMENTS_3=J(4,4,0)

matrix COV_MOMENTS_RAW_1=J(4,4,0)
matrix COV_MOMENTS_RAW_2=J(4,4,0)
matrix COV_MOMENTS_RAW_3=J(4,4,0)

matrix COV_MOMENTS_ANALYTIC_1=J(4,4,0)
matrix COV_MOMENTS_ANALYTIC_2=J(4,4,0)
matrix COV_MOMENTS_ANALYTIC_3=J(4,4,0)

matrix COV_MOMENTS_RAW_ANALYTIC_1=J(4,4,0)
matrix COV_MOMENTS_RAW_ANALYTIC_2=J(4,4,0)
matrix COV_MOMENTS_RAW_ANALYTIC_3=J(4,4,0)

tempfile temp1 temp2 temp3



use "allest", clear

*******************************************
****1) transform data (logs, diff, lag)****
*******************************************

xtset HID YEAR
foreach var of varlist AEQCONS AEQINC {
gen ln`var'=ln(`var')
gen ln`var'lag=l.ln`var'
gen diffln`var'=d.ln`var'
}

gen CONS=AEQCONS*MAEQSIZE
gen INC=AEQINC*MAEQSIZE

egen ln_vill_inc=mean(lnAEQINC), by(VILLAGE YEAR)
egen ln_vill_cons=mean(lnAEQCONS), by(VILLAGE YEAR)
egen vill_inc=mean(AEQINC), by(VILLAGE YEAR)


********************************************
****2) reproduce sum stats *********
********************************************

*A.2, Table 9
*Aurepalle = 1, Kanzara= 2, Shirapur = 3  
foreach i of numlist 1 2 3{
di "now summarizing Village " `i'
su CONS if VILLAGE==`i'
su AEQCONS if VILLAGE==`i'
su INC if VILLAGE==`i'
su AEQINC if VILLAGE==`i'
su MAEQSIZE if VILLAGE==`i'
}


********************************************
****3) reproduce regression ********
********************************************

xtset HID YEAR

***beta_dcdy Section 3, Table 1, Panel B

xi: regress d.lnAEQCONS  d.lnAEQINC i.YEAR if VILLAGE==1 , cluster(HID)
	mat	MOMENTS[2,1]=_b[D.lnAEQINC]
	mat AUX=e(V)
	mat COV_MOMENTS_ANALYTIC_1[2,2]=AUX[1,1]
				est2vec reducedform0, replace vars(d.lnAEQINC) name(Aurepalle)
xi: regress d.lnAEQCONS  d.lnAEQINC i.YEAR if VILLAGE==2 , cluster(HID)
	mat	MOMENTS[2,2]=_b[D.lnAEQINC]
	mat AUX=e(V)
	mat COV_MOMENTS_ANALYTIC_2[2,2]=AUX[1,1]
				est2vec reducedform0, addto(reducedform0) name(Kanzara)
xi: regress d.lnAEQCONS  d.lnAEQINC i.YEAR if VILLAGE==3 , cluster(HID)
	mat	MOMENTS[2,3]=_b[D.lnAEQINC]
	mat AUX=e(V)
	mat COV_MOMENTS_ANALYTIC_3[2,2]=AUX[1,1]
				est2vec reducedform0, addto(reducedform0) name(Shirapur)
					
				
				est2tex reducedform0, preserve path("$latex") levels (90 95 99) fancy replace mark(stars) label

****beta_dcdy without time-demeaning, Section A.11, Table 18
forvalues i=1(1)3{
xi: regress d.lnAEQCONS  d.lnAEQINC  if VILLAGE==`i' , cluster(HID)
	mat	MOMENTSRAW[2,`i']=_b[D.lnAEQINC]
	mat AUX=e(V)
	mat COV_MOMENTS_RAW_ANALYTIC_`i'[2,2]=AUX[1,1]
}
				

***Section 3, Figure 2, Scatter plot of consumption and income growth (after time-demeaning)
**generate residuals for the conditional graph
qui xi: regress d.lnAEQCONS  i.YEAR, cluster(HID)
predict rescons, residuals
qui xi: regress d.lnAEQINC  i.YEAR, cluster(HID)
predict resinc, residuals
twoway (scatter rescons resinc, sort), ytitle(Consumption growth) xtitle(Income growth) title("Conditional on aggregate income", size(medium))	
graph save Graph "$dofiles\Conditional0.gph", replace 
gr export "$latex\icrisat_scatter0_V4.pdf", replace
	
********************************************
****3) estimate the main asymmetry moment***
********************************************
				
***beta_dcdy>0-beta_dcdy<=0 , Section 3, Table 2, Panel B

foreach var of varlist AEQCONS AEQINC {
egen diffln`var'mean=mean(diffln`var'), by(VILLAGE YEAR)
gen devdiffln`var'=diffln`var'-diffln`var'mean
}

gen win=(devdifflnAEQINC>0) if devdifflnAEQINC!=.
gen iawin=win*devdifflnAEQINC

xi: reg devdifflnAEQCONS devdifflnAEQINC win iawin  if VILLAGE==1, cluster(HID)
		est2vec reducedform_winlose0, replace vars(devdifflnAEQINC iawin) name(Aurepalle)
		nlcom (_b[iawin]), post
		mat MOMENTS[4,1]=e(b)
		mat COV_MOMENTS_ANALYTIC_1[4,4]=e(V)
		
		
xi: reg devdifflnAEQCONS devdifflnAEQINC win iawin  if VILLAGE==2	, cluster(HID)
		est2vec reducedform_winlose`kk', addto(reducedform_winlose0) name(Kanzara)
		nlcom (_b[iawin]), post
		mat MOMENTS[4,2]=e(b)
		mat COV_MOMENTS_ANALYTIC_2[4,4]=e(V)
		
		
xi: reg devdifflnAEQCONS devdifflnAEQINC win iawin  if VILLAGE==3	, cluster(HID)
		est2vec reducedform_winlose0, addto(reducedform_winlose0) name(Shirapur)
		nlcom (_b[iawin]), post
		mat MOMENTS[4,3]=e(b)
		mat COV_MOMENTS_ANALYTIC_3[4,4]=e(V)
		
			
		est2tex reducedform_winlose0, preserve path("$latex") levels (90 95 99) fancy replace mark(stars) label
				
			
****beta_dcdy>0-beta_dcdy<=0 , without time-demeaning, Section A.11, Table 18
replace win=d.lnAEQINC>0 if d.lnAEQINC!=.
replace iawin=win*d.lnAEQINC



forvalues i=1(1)3{

xi: reg d.lnAEQCONS d.lnAEQINC win iawin  if VILLAGE==`i', cluster(HID)
		nlcom (_b[iawin]), post
		mat MOMENTSRAW[4,`i']=e(b)
		mat COV_MOMENTS_RAW_ANALYTIC_`i'[4,4]=e(V)				
}				


*************************************************************************
****4) estimate variance and relative variance (and bootstrap)***********
*************************************************************************

***conditional variance of c over conditional variance of y and conditional variance of c for winners and losers 
***Section 3, Table 1, Panel A
forvalues i=1(1)3{
qui su devdifflnAEQCONS if VILLAGE==`i', detail
mat AUX1=r(Var)


qui su devdifflnAEQINC if VILLAGE==`i', detail
mat AUX2=r(Var)

mat MOMENTS[1,`i']=inv(AUX2)*AUX1


qui su devdifflnAEQCONS if VILLAGE==`i' & devdifflnAEQINC>0 & devdifflnAEQINC!=., detail
mat AUX3=r(Var)
qui su devdifflnAEQCONS if VILLAGE==`i' & devdifflnAEQINC<=0 & devdifflnAEQINC!=. , detail
mat AUX4=r(Var)
mat MOMENTS[3,`i']=inv(AUX2)*(AUX3-AUX4)
***Online Appendix, Section A.11, Table 18
qui su d.lnAEQCONS if VILLAGE==`i', detail
mat AUX1=r(Var)


qui su d.lnAEQINC if VILLAGE==`i', detail
mat AUX2=r(Var)

mat MOMENTSRAW[1,`i']=inv(AUX2)*AUX1

qui su d.lnAEQCONS if VILLAGE==`i' & d.lnAEQINC>0 & d.lnAEQINC!=., detail
mat AUX3=r(Var)
qui su d.lnAEQCONS if VILLAGE==`i' & d.lnAEQINC<=0 & d.lnAEQINC!=. , detail
mat AUX4=r(Var)

mat MOMENTSRAW[3,`i']=inv(AUX2)*(AUX3-AUX4)

}


*************************************
***define the AR1MOMENTS ************
*************************************

matrix define lcvcv=J(3,3,0)
matrix define lcvcv_fe_uncond=J(3,3,0)
matrix define lcvcv_fe_cond=J(3,3,0)

*Income process without time-demeaning, Online Appendix, A.11, Table 17, rho and var_alphai
forvalues i = 1(1)3 {
xi: xtreg lnAEQINC lnAEQINClag if VILLAGE==`i', i(HID) fe cluster(HID)
gen sample=e(sample)
***and the variance of the alpha_i
predict u_hat, u
su u_hat if YEAR==2 & VILLAGE==`i', detail
drop u_hat
qui xi: reg lnAEQINC i.HID if VILLAGE==`i' & sample==1
predict reslhs if sample==1, res
qui xi: reg lnAEQINClag  i.HID if VILLAGE==`i' & sample==1
predict resrhs if sample==1, res
reg reslhs resrhs
su resrhs
mat lcvcv_fe_uncond[1,`i']=_b[resrhs]*r(Var)
drop reslhs resrhs sample

}


*Income process benchmark (Household FE and time-demeaned), Section 4.1.2, Table 4, rho and var_alphai
/*Note that the rho we use in the paper is not the regression rho, because we calculate
*rho `by hand' as Cov_y/Var_y, where the variance in the denominator is the variance over 6 rounds. In the regression
*rho, one uses the variance over 5 seasons. The difference is very minor.*/ 
forvalues i = 1(1)3 {
xi: xtreg lnAEQINC lnAEQINClag i.YEAR if VILLAGE==`i', i(HID) fe cluster(HID)
gen sample=e(sample)
**and the variance of the alpha_i
predict u_hat, u
su u_hat if YEAR==2 & VILLAGE==`i', detail
drop u_hat
qui xi: reg lnAEQINC i.HID i.YEAR if VILLAGE==`i' & sample==1
predict reslhs if sample==1, res
qui xi: reg lnAEQINClag  i.HID i.YEAR if VILLAGE==`i' & sample==1
predict resrhs if sample==1, res
reg reslhs resrhs
su resrhs
mat lcvcv_fe_cond[1,`i']=_b[resrhs]*r(Var)
drop reslhs resrhs sample

}


* A.13, Table 21, income process without household fixed effects, just time-demeaned 
foreach var of varlist AEQCONS AEQINC {
egen ln`var'mean=mean(diffln`var'), by(VILLAGE YEAR)
gen devln`var'=ln`var'-ln`var'mean
gen devln`var'lag=l.devln`var'
}

forvalues i = 1(1)3{
corr devlnAEQINC devlnAEQINClag if VILLAGE==`i', covariance 
matrix  lcvcv[1,`i']=r(cov_12)

su d.lnAEQINC if VILLAGE==`i', detail
matrix  lcvcv[2,`i']=r(Var)

su devlnAEQINC if VILLAGE==`i', detail
matrix  lcvcv[3,`i']=r(Var)
}

drop win iawin



save allest0, replace

***this calculates the variance after taking out household fixed effects


*****levels
*Section A.11, Table 17, Var_epsilon, and now you can also calculate rho.
forvalues i=1(1)3 {
xi: reg lnAEQINC i.HID if VILLAGE==`i'

predict res1, res
gen sample1=e(sample)
su res1 if sample1==1, detail
matrix lcvcv_fe_uncond[3,`i']=r(Var)
drop res1 sample1
}

*Section 4.1.2, Table 3, Var_epsilon, and now you can also calculate rho 
forvalues i=1(1)3 {
xi: reg lnAEQINC i.HID i.YEAR if VILLAGE==`i'

predict res1, res
gen sample1=e(sample)
su res1 if sample1==1, detail
matrix lcvcv_fe_cond[3,`i']=r(Var)
drop res1 sample1
}



******and bootstrap to get the standard errors of the variance and ratio of variance

matrix MOMENTS1VARBS=J(1000,4,0)
matrix MOMENTS2VARBS=J(1000,4,0)
matrix MOMENTS3VARBS=J(1000,4,0)


matrix MOMENTS1RAWVARBS=J(1000,4,0)
matrix MOMENTS2RAWVARBS=J(1000,4,0)
matrix MOMENTS3RAWVARBS=J(1000,4,0)

matrix ARMOMENTS1BS=J(1000,3,0)
matrix ARMOMENTS2BS=J(1000,3,0)
matrix ARMOMENTS3BS=J(1000,3,0)

set seed 3
forvalues k=1(1)1000{
di `k'
forvalues i=1(1)3{
use allest0.dta, clear
xtset HID YEAR

keep if VILLAGE==`i'

bsample, cluster(HID) idcluster(NEWHID)
xtset NEWHID YEAR

qui su devdifflnAEQCONS , detail
mat AUX1=r(Var)

qui su devdifflnAEQINC, detail
mat AUX2=r(Var)
*Section 3, Table 1, variance, this is needed for the variance covariance matrix of the moments
mat MOMENTS`i'VARBS[`k',1]=inv(AUX2)*AUX1
gen win=(devdifflnAEQINC>0) if devdifflnAEQINC!=.
gen iawin=win*devdifflnAEQINC

qui su devdifflnAEQCONS  if win>0 & devdifflnAEQINC!=., detail
mat AUX3=r(Var)
qui su devdifflnAEQCONS  if win<=0 & devdifflnAEQINC!=. , detail
mat AUX4=r(Var)
*Section 3, variance difference
mat MOMENTS`i'VARBS[`k',3]=inv(AUX2)*(AUX3-AUX4)

****without time-demeaning
qui su d.lnAEQCONS , detail
mat AUX1=r(Var)

qui su d.lnAEQINC, detail
mat AUX2=r(Var)
*Section A.11, Table 18, standard errors for ratio of variance
mat MOMENTS`i'RAWVARBS[`k',1]=inv(AUX2)*AUX1

qui su d.lnAEQCONS  if d.lnAEQINC>0 & d.lnAEQINC!=., detail
mat AUX3=r(Var)
qui su d.lnAEQCONS  if d.lnAEQINC<=0 & d.lnAEQINC!=. , detail
mat AUX4=r(Var)
*Section A.11, Table 18, bootstrapping standard errors for difference of variance
mat MOMENTS`i'RAWVARBS[`k',3]=inv(AUX2)*(AUX3-AUX4)


****REGRESSION MOMENTS
*Section 3, Table 1 and 2, this is to calculate the moment covariance matrix
qui xi: regress d.lnAEQCONS  d.lnAEQINC i.YEAR, cluster(HID)
	mat	MOMENTS`i'VARBS[`k',2]=_b[D.lnAEQINC]
 
qui xi: reg devdifflnAEQCONS devdifflnAEQINC win iawin , cluster(HID)
	nlcom (_b[iawin]), post
	mat MOMENTS`i'VARBS[`k',4]=e(b)
}
}



***AR MOMENTS BOOTSTRAPPING


forvalues k=1(1)1000{
di `k'
forvalues i=1(1)3{
use allest0.dta, clear
xtset HID YEAR

keep if VILLAGE==`i'

bsample, cluster(HID) idcluster(NEWHID)
xtset NEWHID YEAR

xi: xtreg lnAEQINC lnAEQINClag i.YEAR if VILLAGE==`i', i(NEWHID) fe cluster(HID)
gen sample=e(sample)
***and the standard error of the alpha_i
mat ARMOMENTS`i'BS[`k',2]= e(sigma_u)

***and for the variance of rho
qui xi: reg lnAEQINC i.HID i.YEAR if VILLAGE==`i' & sample==1
predict reslhs if sample==1, res
qui xi: reg lnAEQINClag  i.HID i.YEAR if VILLAGE==`i' & sample==1
predict resrhs if sample==1, res
reg reslhs resrhs
su resrhs
matrix ARMOMENTS`i'BS[`k',1]=_b[resrhs]*r(Var)
drop reslhs resrhs sample

**and for the variance of the income. 
xi: reg lnAEQINC i.NEWHID i.YEAR if VILLAGE==`i'
predict res1, res
gen sample1=e(sample)
su res1 if sample1==1, detail
matrix ARMOMENTS`i'BS[`k',3]=r(Var) 
drop res1 sample1

	
	
}
} 

svmat MOMENTS1VARBS, names(moments1varbs)
svmat MOMENTS2VARBS, names(moments2varbs)
svmat MOMENTS3VARBS, names(moments3varbs)

svmat MOMENTS1RAWVARBS, names(moments1rawvarbs)
svmat MOMENTS2RAWVARBS, names(moments2rawvarbs)
svmat MOMENTS3RAWVARBS, names(moments3rawvarbs)

svmat ARMOMENTS1BS, names(armoments1bs)
svmat ARMOMENTS2BS, names(armoments2bs)
svmat ARMOMENTS3BS, names(armoments3bs)

svmat ARMOMENTS1BS, names(armoments_table_1bs)
svmat ARMOMENTS2BS, names(armoments_table_2bs)
svmat ARMOMENTS3BS, names(armoments_table_3bs)
*turn the covariance into rho
replace armoments_table_1bs1=armoments1bs1/armoments1bs3
replace armoments_table_2bs1=armoments2bs1/armoments2bs3
replace armoments_table_3bs1=armoments3bs1/armoments3bs3
*turn the standard deviation of the fixed effects into a variance
replace armoments_table_1bs2=armoments1bs2^2
replace armoments_table_2bs2=armoments2bs2^2
replace armoments_table_3bs2=armoments3bs2^2
***and turn the variance of income into the variance of errors
replace armoments_table_1bs3=armoments1bs3*(1-armoments1bs1^2)
replace armoments_table_2bs3=armoments2bs3*(1-armoments2bs1^2)
replace armoments_table_3bs3=armoments3bs3*(1-armoments3bs1^2)

su armoments_table_1bs*
su armoments_table_2bs*
su armoments_table_3bs*

forvalues i=1(1)3 {
su moments`i'varbs1
mat COV_MOMENTS_`i'[1,1] = r(sd)^2
su moments`i'varbs3
mat COV_MOMENTS_`i'[3,3] = r(sd)^2
su moments`i'varbs2
mat COV_MOMENTS_`i'[2,2] = r(sd)^2
su moments`i'varbs4
mat COV_MOMENTS_`i'[4,4] = r(sd)^2

su moments`i'rawvarbs1
mat COV_MOMENTS_RAW_`i'[1,1] = r(sd)^2
su moments`i'rawvarbs3
mat COV_MOMENTS_RAW_`i'[3,3] = r(sd)^2
su moments`i'rawvarbs2
mat COV_MOMENTS_RAW_`i'[2,2] = r(sd)^2
su moments`i'rawvarbs4
mat COV_MOMENTS_RAW_`i'[4,4] = r(sd)^2

forvalues k1=1(1)4 {
	forvalues k2=1(1)4 {
	if `k1'!=`k2' {
	corr moments`i'varbs`k1' moments`i'varbs`k2', cov 
	mat COV_MOMENTS_`i'[`k1',`k2'] = r(cov_12)
	
	corr moments`i'rawvarbs`k1' moments`i'rawvarbs`k2', cov 
	mat COV_MOMENTS_RAW_`i'[`k1',`k2'] = r(cov_12)
	}
	}
	} 


} 


drop _all
svmat MOMENTS, names(moments)
svmat MOMENTSRAW, names(momentsraw)

svmat COV_MOMENTS_1, names(cov_moments_1)
svmat COV_MOMENTS_2, names(cov_moments_2)
svmat COV_MOMENTS_3, names(cov_moments_3)

svmat COV_MOMENTS_RAW_1, names(cov_moments_raw_1)
svmat COV_MOMENTS_RAW_2, names(cov_moments_raw_2)
svmat COV_MOMENTS_RAW_3, names(cov_moments_raw_3)

svmat COV_MOMENTS_ANALYTIC_1, names(cov_moments_analytic_1)
svmat COV_MOMENTS_ANALYTIC_2, names(cov_moments_analytic_2)
svmat COV_MOMENTS_ANALYTIC_3, names(cov_moments_analytic_3)

svmat COV_MOMENTS_RAW_ANALYTIC_1, names(cov_moments_raw_analytic_1)
svmat COV_MOMENTS_RAW_ANALYTIC_2, names(cov_moments_raw_analytic_2)
svmat COV_MOMENTS_RAW_ANALYTIC_3, names(cov_moments_raw_analytic_3)



order moments1* moments2* moments3* moments4* cov_moments_1* cov_moments_2* cov_moments_3* cov_moments_analytic* momentsraw* cov_moments_raw*
outsheet moments1* moments2* moments3* moments4* cov_moments_1* cov_moments_2* cov_moments_3* cov_moments_analytic* momentsraw* cov_moments_raw*  using "$output/icrisatmoments0", nonames replace


drop _all
svmat lcvcv, names(armoments)
outsheet armoments* using "$programmes/armoments", nonames replace
outsheet armoments* using "$output/armoments", nonames replace

svmat lcvcv_fe_cond, names(armoments_fe_cond)
outsheet armoments_fe_cond* using "$programmes/armoments_fe0", nonames replace
outsheet armoments_fe_cond* using "$output/armoments_fe0", nonames replace

svmat lcvcv_fe_uncond, names(armoments_fe_uncond)
outsheet armoments_fe_uncond* using "$programmes/armoments_fe1", nonames replace
outsheet armoments_fe_uncond* using "$output/armoments_fe1", nonames replace

*****and the village-specific data set for the bivariate correlations
use allest0.dta, clear
preserve
sort HID YEAR 
keep if VILLAGE==2
keep if lnAEQCONS!=. 
xi: xtreg lnAEQCONS i.YEAR if VILLAGE==2, i(HID) fe 
predict reslnAEQCONS if e(sample)==1, res
sort HID YEAR
outsheet HID YEAR lnAEQCONS reslnAEQCONS using "$output/perm_vill_2" if VILLAGE==2, nonames replace
restore

